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ABSTRACT 

Numerous studies have investigated the role of thermal instability in regulating the phase transition between 
the cold cloudy and warm diffuse medium of the interstellar medium. Considerable interest has also been 
devoted to investigating the properties of turbulence in thermally unstable flows, with a special emphasis on 
molecular clouds and the possibility of star formation. In this study, we investigate another setting in which this 
instability may be important, namely its effect on dynamo action in interstellar flows. The setup we consider is 
a three dimensional periodic cube of gas with an initially weak magnetic field, subject to heating and cooling, 
the properties of which are such that thermal instability is provoked in a certain temperature regime. Dynamo 
action is established through external forcing on the flow field. By comparing the results with a cooling function 
with exactly the same net effect but no thermally unstable regime, we find the following. Reference runs with 
non-helical forcing were observed to produce no small-scale dynamo action below the Reynolds number 97. 
Therefore, we expect the magnetic fields generated in the helical runs to be purely due to the action of a 
large-scale dynamo mechanism. The critical Reynolds number for the onset of the large-scale dynamo was 
observed to roughly double between the thermally stable versus unstable runs, the conclusion being that the 
thermal instability makes large-scale dynamo action more difficult. Whereas density and magnetic fields were 
observed to be almost completely uncorrected in the thermally stable cases investigated, the action of thermal 
instability was observed to produce a positive correlation of the form B oc p 2 . This correlation is rather weak, 
and in addition it was observed to break down at the limit of highest densities. 

Subject headings: magnetohydrodynamics, instabilities, turbulence, ISM: general 



1. INTRODUCTION 

The role of thermal instability (hereafter TI) in regulating 
the phase transition between the cold and warm neutral com- 
ponents, normally denoted with CNM and WNM, of the in- 
terstellar matter (hereafter ISM) has been intensively studied 
over several decades (for a review, see e.g. Cox 2005) since 
the original seminal paper by Field, Goldsmith, & Habing 
(1969) proposed this mechanism, known as the FGH-model. 
According to it, two thermally stable phases (cold and cloudy; 
warm and diffuse) co-exist in pressure equilibrium regulated 
by the presence of a thermally unstable phase at an interme- 
diate temperature. Now it is observationally well-established 
that the picture is not quite simple as that, but there is also a 
third, hot and dilute, component of the ISM. The first model 
that incorporated the hot component, generated by the action 
of supernova explosions (hereafter SNe), to the FGH model 
was presented by McKee & Ostriker (1977). This model sug- 
gested rather drastic modifications to the structure of ISM: 
with the inclusion of SN-action, most (70-80 percent) of the 
gas became filled by the hot component. Since then, with 
extensive numerical modelling of the SN-driven flows (e.g. 
Chiang & Prendergast 1985; Rosen et al. 1993; Korpi et al. 
1999; Gazol-Patino & Passot 1999; de Avillez 2000; Balsara 
et al. 2004; Gressel et al. 2008) and more careful observa- 
tional determinations (for a review, see e.g. Ferriere 2001), 
the estimates of the filling factor of the hot component have 
been reduced to 10-30 percent near the Galactic midplane, 
being bigger at larger heights. Moreover, most of the hot gas 
seems to be confined in large bubbles created by clustered SN 
activity rather than being distributed homogeneously around 
the Galaxy (e.g. Ferriere 2001). 

A number of investigations of interstellar turbulence at rel- 
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evant scales for studying molecular cloud formation or any 
other process of interest in the atomic ISM still use the FGH- 
picture as a starting point. This seems justified in the light of 
the hot gas filling factor being relatively small near the Galac- 
tic midplane. Various numerical studies investigating the in- 
teraction of turbulence and TI have been published. In most of 
these studies turbulence is forced by sources other than the TI 
itself: random turbulent forcing at varying scales and Mach 
numbers (e.g. Gazol et al. 2005), localized injections of en- 
ergy mimicking stellar winds (e.g. Vazquez-Semadeni et al. 
2000), the magnetorotational instability (e.g. Piontek & Os- 
triker 2007), and systematic large-scale motions such as prop- 
agating shock fronts (e.g. Koyama & Inutsuka 2002) and con- 
verging flows e.g. by Audit & Hennebelle (2005), Heitsch 
et al. (2005) and Vazquez-Semadeni et al. (2006) have been 
considered. One of the major findings from these models is 
that, due to the turbulence present in the system, large pres- 
sure deviations are generated and significant amounts of gas 
can exist in the thermally unstable regime. These results sug- 
gest, as already noted by Norman & Ferrara (1996), that the 
FGH picture of the ISM exhibiting 'discrete' temperatures and 
densities and a unique equilibrium pressure should be modi- 
fied in the direction of a 'continuum' of states with an overall 
pressure balance but with large deviations from it. 

In recent years the possibility of driving turbulence by the 
TI itself has also received some attention. Contrary to Krit- 
suk & Norman (2002), who found turbulence to die out as a 
power law, Inutsuka & Koyama (2007) found the turbulence 
to be sustained — at least for times up to 0.1 Gyr, although the 
turbulent motions were too weak to be significant in the ISM. 
However, in the study by Brandenburg, Korpi, & Mee (2007), 
hereafter BKM07, the Tl-generated turbulence was observed 
to decay in simulations extending to more than 1 Gyr. Even 
the inclusion of large-scale shear did not help sustain the tur- 
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bulence. Therefore, it seems evident that TI alone cannot 
produce self-sustained turbulence, in contrast to the gravita- 
tional instability studied by Gammie (2001) or the magne- 
torotational turbulence studied e.g. by Hawley, Gammie, & 
Balbus (1995) and Brandenburg et al. (1995). 

The line-of-sight magnetic field strengths in dense molecu- 
lar clouds, as revealed by the Zeeman splitting measurements 
of molecular lines, follow an approximate correlation of the 
form B|| oc n(H 2 ) 1 ^ 2 , where n(H 2 ) is the density of the 
molecular hydrogen in the cloud, e.g. Crutcher (1999) (for 
a recent review, see e.g. Heiles & Crutcher 2005). For dif- 
fuse clouds at lower densities n ~ 0.1 - 100 cm 3 also seen 
through the Zeeman splitting of the HI line, no clear correla- 
tion between the magnetic field strength and density has been 
found (Troland & Heiles 1986; Myers, Goodman, Gusten, & 
Heiles 1995; Heiles & Crutcher 2005). Passot & Vazquez- 
Semadeni (2003) suggested that this phenomenon could be 
due to the strongly varying Alfvenic Mach number, Ma a , in 
the turbulent ISM. At low values of Ma a , their analysis re- 
vealed an anti-correlation between the density and magnetic 
pressure due to the dominance of slow MHD waves. At higher 
Ma a , both slow and fast modes were possible, their interfer- 
ence leading to de-correlation for intermediate densities, and 
re-correlation for the highest densities. Heitsch et al. (2004) 
proposed an alternative mechanism, namely ambipolar dif- 
fusion, enhanced by the vigorous interstellar turbulence to 
be rapid enough to de-correlate the magnetic field from the 
density. The fast reconnection of magnetic field lines in the 
cloudy medium has also been discussed in this context, e.g. 
by Santos-Lima et al. (2010). 

In addition, there is observational evidence of it being par- 
ticularly the regular component, seen through synchrotron 
emission, that is de-correlated from the gas density. For 
instance, observations of barred spiral galaxies indicate the 
magnetic fields forming coherent structures over the length 
scale of the bars themselves are not reacting as strongly as 
expected to the compression by the systematic motions oc- 
curring in the bars (e.g. NGC 1097 and NGC 1365 analyzed 
by Beck et al., 2005). A similar situation has been reported 
in M51 by Fletcher et al. (2010), where the regular magnetic 
field is not strongest at the locations of the highest compres- 
sion due to the spiral density wave. It has been proposed that 
such a loose or completely nonexistent correlation between 
the regular component of the magnetic field and the diffuse 
clouds could be a result of the detachment of the large-scale 
field from the clouds in a fragmentation process such as the 
one resulting from thermal or gravitational instability (Beck 
et al. 2005). 

During the recent years, numerical models in the magne- 
tohydrodynamic regime including thermal and gravitational 
instabilities able to reach down to the size of the molecular 
clouds themselves have been developed by e.g. Hennebelle 
et al. (2008), Heitsch & Hartmann (2008) and Banerjee et al. 
(2009). The results from these calculations seem to be in very 
good agreement with the observations, i.e. producing very 
weak correlation between the magnetic field and density for 
lower densities, while re-establishing the correlation of the 
form B oc n 1 / 2 at higher densities. Our first motivation in 
this manuscript is to isolate the role of thermal instability in 
the process that produces the density-magnetic field correla- 
tion. 

Another interesting question that has not gained much at- 
tention is the effect of thermal instability on the galactic dy- 



namo mechanism. Evidently, thermal instability cannot be 
expected by itself to drive turbulence that could generate and 
maintain a self-sustained dynamo. Rather, it can be expected 
to influence the dynamics of a flow that is dynamo active, and 
therefore may also have an effect on the properties of the dy- 
namo. This is the second aspect we are aiming at investigating 
in detail in the present manuscript. 

We adopt a very simple approach: as in our previous study 
of thermal instability under a purely hydrodynamic setting 
(BKM07), we use a cubic periodic computational domain 
filled with an initially uniform, non-stratified gas that has 
the characteristics of the atomic component of the interstel- 
lar matter. Unlike BKM07, a weak random magnetic field is 
inserted into the system at the initial state. Dynamo action 
is achieved by using an external forcing function on the flow 
field. In the case where the forcing is helical, a large-scale 
dynamo develops. For non-helical forcings, only small-scale 
dynamo action can be expected in our setting. The system is 
then subjected to either thermally stable or unstable cooling 
functions balanced by a constant UV-heating function. From 
the resulting set of models, we monitor the growth rate and 
saturation properties of the different types of dynamos, and 
the correlation between the density and magnetic field devel- 
oping in each system. 

2. MODEL 

2.1. Governing equations 

We consider the governing equations for a compressible 
and magnetized perfect gas, 
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where u is the velocity, p the density, s the specific entropy, 
A the vector potential of the magnetic field B = V x A, J = 
^o^V x B the current, with Sy = ^{ui,j+Uj,i) — I^V-m 
being the traceless rate of strain tensor, v is the kinematic vis- 
cosity, x is the thermal diffusivity, rj the magnetic diffusiv- 
ity, /Uo the permeability of free space, and C is the net cool- 
ing/heating, i.e. the difference between cooling and heating 
functions, with 

C = P A-T, (5) 

where V = const is assumed for the heating function; here 
we consider the photoelectric heating by interstellar grains 
caused by the stellar UV radiation field, for which Wolfire et 
al. (1995) give the value of 0.015 ergg -1 s _1 at n = 1 cm~ 3 . 

Following common practice, we adopt a perfect gas where 
p and s are related to pressure p and temperature T via the 
relations 



P 



K nn 

-PT, 
I 1 



s = c v Inp — c p \np + so, 



(6) 



where 1Z = 8.314 x 10 7 cm 2 s _2 K _1 is the universal gas 
constant, and p is the mean molecular weight; in this study 
we have adopted p = 0.62 in all cases, corresponding to 
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complete ionization. This is naturally unrealistic for the cold 
and cool phases, for which the ionization fraction can be ex- 
pected to be low. To choose a higher mean molecular weight 
better corresponding to the cooler phases, however, would 
be equally unrealistic, as a significant fraction of the gas is 
still expected to be ionized. A better approach to deal with 
this issue would be to follow the ionization fraction and ad- 
just n accordingly. This, however, is out of the scope of the 
present study. The rest of the thermodynamic quantities read 
1Z//J, = c p — c v , with c p and c v being the specific heats at 
constant pressure and volume, respectively, 7 = c p /c v = 5/3 
is their assumed ratio, c s is the adiabatic sound speed, T is the 
temperature, and the two are related to the other quantities via 
c 2 = jlZTf/j,. The specific entropy is defined up to a constant 
so, the value of which is unimportant for the dynamics. 

We adopt a parameterization of the cooling function equal 
to that given by Sanchez-Salcedo et al. (2002), which has been 
obtained by fitting a piecewise power law function of the form 

A(T) = C M+1 T^-+ 1 for Ti <T <T i+u (7) 

to the equilibrium pressure curve of the standard model of 
Wolfire et al. (1995) for the ISM in the solar neighborhood. 
This cooling function has a thermally unstable branch in be- 
tween the temperatures T = 313 K and 6102 K, where /3 < 1; 
henceforth, we denote this cooling function with 'SS'. To be 
able to extract the effect of the thermally unstable branch of 
the cooling function on the flow dynamics, we fabricate a ther- 
mally stable counterpart, denoted with 'TS', with the criterion 
that for constant density, the net cooling caused by these two 
functions is the same. The coefficients for both cooling func- 
tions are given in Table 1, and a comparison plot is presented 
in Fig. 1. In this figure, we plot the heating-cooling equilib- 
rium pressures with the two cooling functions, with respect to 
the chosen value of the UV-heating function T. 

We consider the case where turbulence is driven by an addi- 
tional body force, /force, in the momentum equation. Follow- 
ing Kapyla & Brandenburg (2009), we use a forcing function 
/force given by 

/forced, t) = Rc{Nf k(t) exp[tk(t) ■ x + i<j>(t)}} , (8) 

where x is the position vector, N — foc s (kc s / St) 1 / 2 is a 
normalization factor, /o is the forcing amplitude, k = |fc| 
the wavevector, St is the length of the time step, and — 7r < 
4>(t) < tt a random delta-correlated phase. The wavevector 
and the random phase change at every time step, resulting in 
a forcing function that is J-correlated in time. We force the 
system with transverse helical waves, 

f k = r ■ /r hcl with r 13 = 5 ^-_^ k \ (9 ) 

V 1 + <? 

The helical case with positive helicity corresponds to a = 1, 
with the forcing function 

/r hcl = (k x e) /^/k 2 -(k-ef, (10) 

where e is an arbitrary unit vector not aligned with the forcing 
wavevector. 

Thus, //, can be chosen to give either non-helical (a = 0) 
or helical (cr = 1) transversal waves with \fk\ 2 = 1, where k 
is chosen randomly from a predefined range in the vicinity of 
the average non-dimensional wavenumber kf / k\ at each time 
step. Here k\ is the wavenumber corresponding to the domain 



TABLE 1 

Coefficients for the thermally unstable (SS) and stable (TS) 
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size, and k[ is the wavenumber of the energy-carrying scale, 
chosen to be between 2.5 and 3.5 times the smallest wavenum- 
ber in the box, fci = 27r/(0.2 kpc). In a system very similar to 
the one under study here, namely Brandenburg (2001), large- 
scale dynamo action was found with helical forcing. In the 
presence of shear, non-helical forcing has also been found to 
excite a large-scale dynamo (e.g. Yousef et al. 2008), but in 
the absence of it, can lead only to small-scale dynamo action 
(e.g. Haugen et al. 2004). In this study, both types of forcings 
are investigated, in the hope to study the effect of TI both on 
the small- and large-scale dynamos. We systematically vary 
the type and vigor of the forcing, and subject the gas either to 
thermally stable or unstable cooling function; the parameters 
for the runs produced can be found from Table 2. 

It is convenient to measure time in gigayears (Gyr), speed in 
km/s, and density in units of 10~ 24 g cm~ 3 . The resulting unit 
of pressure is 10~ 14 dyn, and the unit of length 1( km/s) x 
Gyr s=y 1 kpc. Viscosity and thermal diffusivity are measured 
in units of Gyr km 2 s~ 2 . The resulting unit for energy is 
Eo=3.2x 10 50 ergs, and the unit of magnetic field is 0.35/iG. 

We use periodic boundary conditions in all three directions 
for a computational domain of size (200 pc) 3 . Near the galac- 
tic midplane, the cell size reported from SN-turbulence range 
in between 60-100pc (e.g. Korpi et al. 1999); this can be 
thought as the relevant length scale for the galactic dynamo 
near the midplane, one of the basic ingredients of which SN- 
turbulence is thought to be (see e.g. Ferriere 2001). There- 
fore, the chosen computational domain size can be thought of 
as a compromise in being large enough to contain a few turbu- 
lent SN cells but small enough to follow the phase segregation 
process into cold cloudy and warm dilute phases of the ISM. 

We use the Pencil Code, 2 which is a non-conservative, 
high-order, finite-difference code (sixth order in space and 
third order in time) for solving the compressible hydrody- 
namic equations. Because of the non-conservative nature 
of the code, diagnostics giving the total mass and total en- 
ergy (accounting for heating/cooling terms) are monitored 
and simulations are only deemed useful if these quantities are 
in fact conserved to reasonable precision. The mesh spac- 
ings in the three directions are assumed to be the same, i.e. 
Sx = Sy = Sz. In the simulations reported in this manuscript, 
the size of the computational grid is 128 3 , yielding a resolu- 
tion of roughly 1.6pc in our model. We note that a resolution 
study in a corresponding hydrodynamic setup was made in 
our previous study (BKM07), where doubling the resolution 
was observed to not significantly change the flow properties. 

As in our previous study (BKM07), see also Koyama & 
Initsuka (2004) and Piontek & Ostriker (2005), we include 
thermal conduction, which stabilizes the gas at wavelengths 
larger than the critical one of the condensation mode (Field 

2 http://http://code.google.com/p/pencil-code/ 
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1965). This wavelength is usually referred to as the Field 
length, and is three orders of magnitude larger in our sim- 
plified model than in reality (of the order of O.OOlpc). The 
required stabilizing effect could, in principle, be achieved by 
high enough numerical diffusion arising from the scheme it- 
self (see e.g. Gazol et al. 2005), but with our high-order finite- 
difference scheme this effect is not large enough with the cho- 
sen viscosity coefficients. We emphasize that no shock or hy- 
perviscosity has been used in the present simulation. There- 
fore, the only means of stabilizing the code is through the 
regular viscosity v, magnetic diffusivity r], and thermal dif- 
fusivity x- In order to damp unresolved ripples at the mesh 
scale Sx in the trail of structures moving at speed U, the min- 
imum viscosity and minimum diffusion must be of the order 
of about 0.01 x USx (see Brandenburg & Dobler 2002); for 
the setup presented in this paper, the minimal viscosity coeffi- 
cients to resolve TI producing maximal velocities of the order 
of 100 km s _1 with the grid and guarantee numerical stability 
are v, \ and ry of 5 x 10~ 3 , meaning that the Prandtl numbers 
are of unity. The kinetic and magnetic Reynolds numbers are 
defined as 

Re = ii rms Z/ jv and Rm = u ruls lf/r] (11) 

where I / is the length scale corresponding to the external forc- 
ing. Keeping Prandtl numbers as unity, the kinetic and mag- 
netic Reynolds numbers are equal. 

2.2. Diagnostic for the magnetic flux 

Instead of the commonly used point-to-point correlation of 
density and magnetic field strength at each cell of the compu- 
tational grid, we propose to measure the correlation between 
the magnetic field strength and density using a different pro- 
cedure that cleans out the unnecessary information of how 
the data points cluster about the most probable magnetic field 
strength, while revealing more clearly the mean trend in the 
density-magnetic field correlation. First, we divide the loga- 
rithmic density field into n bins of equal size. Then we calcu- 
late the total magnetic flux, J s B ■ dS, through the horizontal 

($tot)> and tne two vert i ca l ($tot and ^"tot) surfaces, and the 
magnetic flux through these surfaces for each density interval 
($^ y , <I>* Z , and ^f 2 ), i = 1, ...,n. To prevent the flux diag- 
nostic being dependent on the amount of gas in each density 
bin, we need to scale with the filling factor of each bin, fa, i.e. 

<J>kl 

m kl = l (12) 

.M^tot 

In principle, both the fluxes and the filling factors are chang- 
ing as functions of time, especially in the dynamo-active runs 
during the growth phase of the magnetic field. Therefore, 
the flux diagnostic is calculated only after the magnetic field 
growth has saturated, after which practically no change in any 
of the quantities can be detected. All the N = {nx, ny, nz} 
surfaces along the three line-of-sights for the corresponding 
{\[>2/ z ) \fy xz j vj/^J*} are added up in the calculation. Comparing 
the fluxes calculated over all the three possible directions also 
allows us to detect anisotropics; this is, however, not a rele- 
vant issue for the present study excluding stratification, rota- 
tion and shear. This diagnostic is expected to show a mono- 
tonically increasing trend for any system showing correlation 
between the density and magnetic field, monotonically de- 
creasing trend for anti-correlated fields, and flat behavior, the 
value approaching unity, if the quantities remain completely 



uncorrected. Before applying this diagnostic to the dynamo- 
generated fields, we test it with artificially generated magnetic 
fields on top of the density field from a corresponding hydro- 
dynamic simulation (see next section). 

3. RESULTS 

The produced runs, their basic parameters, and some 
ensemble-averaged quantities describing the basic results are 
listed in Table 2. Rms velocity u rms , rms total magnetic field 
-Bjms> thermal E t h, kinetic £kin, magnetic E mm energies, 
rms Mach number Ma rms = u rms /c rms , where c rms is the 
rms sound speed, and magnetic Reynolds number Rm, equal 
to Re, are all calculated as averages over the whole volume 
and over several Gyrs of time after the saturated state has been 
reached. For the hydrodynamic runs, instead of the magnetic 
Reynolds number, the fluid Reynolds number is given. The 
growth rate, Ab, is measured from -B*ms during the exponen- 
tial phase of the magnetic field growth; in the table, we list 
its inverse to illustrate the growth time scale of the magnetic 
field in Myrs. To measure the density contrast, Ap, averages 
of the extremum values of densities over the saturated state 
are used. The labelling system of the runs is as follows: 'HD' 
refers to hydrodynamic runs excluding magnetic field; 'TS' 
refer to the usage of thermally stable cooling function; 'SS' 
refers to the thermally unstable cooling function; 'nh' refers 
to the usage of a non-helical forcing function. 

First we demonstrate the operation of the diagnostic de- 
fined in Eq. (12), using some purely hydrodynamic calcula- 
tions otherwise identical to the magnetohydrodynamic ones, 
but excluding the magnetic field. We use the density fields 
from the thermally stable run HDTSa and its thermally unsta- 
ble counterpart HDSSa, on top of which we artificially add a 
magnetic field directed in the positive y-direction. The am- 
plitude of the field is given either by (a) |B| = B p 1/2 , (b) 
\B\ = Bqp^ 1 ! 2 or (c) the amplitude of the field is a random 
number in the range [0, B ], with B fixed to unity. In such 
a setup, only the flux component ^ xz is nonzero, and this 
quantity is plotted for the different cases in Fig. 2. As evident 
from this figure, for a setting with a positive correlation be- 
tween the density and magnetic field strength (case a; plotted 
with solid lines of varying color), the flux diagnostic shows a 
monotonically increasing trend. For an anti-correlated setup 
(case b; plotted with dashed lines of varying color), the diag- 
nostic shows a monotonically decreasing trend. For the case 
of completely uncorrected magnetic field strength and den- 
sity (case c; plotted with dotted lines of varying color) the 
curve is flat and has a value around unity. The density range 
in the thermally unstable Run HDSSa is larger than in the 
stable counterpart, Run HDTSa, but the diagnostic behaves 
similarly in both cases. 

3.1. Hydrodynamic state of the system 

Let us first compare the purely hydrodynamic runs with 
different cooling functions and forcings. For the runs with 
thermally stable cooling functions (labeled with HDTS in Ta- 
ble 2), the formation of density structures occurs solely by 
the mixing and compression caused by the turbulent flow. For 
cases with the weakest forcing, Run HDTSa, rms-velocity of 
the order of 2kms _1 , and density contrast of 3 is observed. 
The density field of this run at 250 Myrs is illustrated in Fig. 3 
leftmost panel, from which it can easily be seen that the tur- 
bulence is so weak that the density structures generated by 
the flow do not show large contrasts and remain quite large in 
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Results from various runs. The thermally stable runs are labelled with 'TS\ thermally unstable with 'SS', purely 
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FIG. 1. — The equilibrium pressures calculated as function of density for 
the unstable 'SS' curve (solid) and stable 'TS' curve (dashed). 

size. When the strength of the forcing is increased, the den- 
sity contrast gets larger, as expected, so that in Run HDTSb, 
with 2.5 times stronger forcing, the contrast is roughly 20, 
and in Run HDTSc, with 3.5 times higher forcing, roughly 
100, the respective rms-velocities being 6 and 8kms _1 . The 
appearance of the flow field, illustrated in Fig. 3, the second 
panel from the left for Run HDTSc, has somewhat changed, 
the density structures generated now showing up as narrow 
filaments. 

The effect of increasing forcing is also apparent in the prob- 
ability distribution functions (hereafter PDFs) of the quanti- 
ties, plotted in Fig. 4 for Run HDTSa (dashed-dotted lines) 
and Run HDTSc (dashed lines). Evidently, the distributions 
of all the quantities are all unimodal, with the width of them 
increasing as the forcing is increased. The density PDF is 
clearly reflecting the somewhat increased capability of the 
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FIG. 2. — Diagnostic Eq. (12) calculated by imposing an artificially gener- 
ated magnetic field on top of the thermally stable hydrodynamic Run HDTSa 
(black lines) and thermally unstable Run HDSSa (orange/grey lines). Differ- 
ent types of magnetic field dependencies on density are plotted with different 
line styles: solid line, |£?| = Bap 1 !' 1 , dashed line. |B| = Bop -1 ! 2 , dotted 
line: random magnetic field. 

flow to generate denser structures with a smaller filling fac- 
tor. The denser the structures the flow is capable of creating, 
the more will they cool, as the cooling term is of the form 
pA, modifying the rather symmetrical temperature PDF for 
Run HDTSa to exhibit an extended wing towards lower tem- 
peratures for Run HDTSc. This excess cooling will aid the 
dense regions in becoming denser, and a similar asymmetrical 
wing towards higher densities is seen in the density PDF for 
Run HDTSc. This is in agreement with the results of Passot 
& Vazquez-Semadeni (1998), who studied highly compress- 
ible polytropic flows with a varying Mach number. For poly- 
tropic indeces less than one, a very similar wing, approaching 







FIG. 3. — Snapshots of the logarithmic density fields from different hydrodynamic runs at 250 Myrs. From left to right: HDTSa, HDTSc, HDSSa, HDSSc. 
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FIG. 4. — Probability density functions for the thermally unstable hydrodynamic runs HDSSc (solid line) and HDSSa (dotted line), and the thermally stable 
hydrodynamic runs HDTSc (dashed line) and HDTSa (dashed-dotted line) at 250 Myrs of the simulation runs. On the left we show the PDF of logarithmic 
density, in the middle the PDF of logarithmic temperature, and on the right the PDF of logarithmic pressure. 
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FIG. 5 . — Snapshots of the magnetic field strength from the thermally stable run TSa and thermally unstable run SSa at 2.5 Gyrs of evolution. From left to right: 
By in TSa, B z in TSa, B y in SSa, and B z in SSa. The minimum value plotted with blue color corresponds to -1.2 fiG and the maximum value to 1.2 fj,G, used 
in all plots. 

to enhance the condensation process, so that excess cooling 
due to turbulent compression leads to an even larger density 
contrast and smaller structures. 

We also plot the PDFs of the quantities for Run HDSSc 
(solid lines) and Run HDSSa (dotted lines) in Fig. 4. For 
the Run HDSSa with relatively weak forcing, the density and 
temperature distributions are still clearly bimodal, as has com- 
monly been observed for thermally unstable flows, e.g. by 
BKM07, although large amounts of gas can be seen in the 
'forbidden' regime of thermally unstable temperatures with 
j3 < 1. For the Run HDSSc, the bimodality of the density and 
temperature PDFs is much less pronounced; both distributions 
have become wider, and the peaks have been shifted towards 
higher densities and lower temperatures. This is reflected also 
in the peak of the pressure distribution. For the strongest forc- 
ing, the most probable pressure is somewhat higher than for 
the weaker forcing. A similar trend in the pressure distribution 
can also be seen in the thermally stable runs. The filling factor 
of the cold component gets considerably smaller as the forc- 



power-law behavior, was seen to develop towards higher den- 
sities. Similar, skewed density PDFs towards higher densities 
have also been reported from 3D-simulations of polytropic 
flows by Nordlund & Padoan (1999), interpreted to be caused 
by the effective 7 being reduced to below one in the regions 
of high compression. 

For the runs with thermally unstable cooling functions (la- 
beled with HDSS in Table 2), a major part of the formation of 
density structures is expected to occur due to the thermal in- 
stability, although the turbulence can also be expected to alter 
this picture considerably (see e.g. BKM07; Gazol et al. 2005). 
As can be seen from the two rightmost panels of Fig. 3, the 
appearance of the flow is very different from the thermally sta- 
ble cases (the two leftmost panels in that figure). The action of 
the thermal instability leads to a very rapid (timescale of few 
tens of Megayears) segregation of the gas into cold cloudy 
and warm diffuse phases, the dense clouds clearly visible in 
the snapshots of the density field. Similarly to the thermally 
stable cases, the increasing vigor of turbulence can be seen 
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ing is increased, especially clearly seen when comparing the 
two leftmost panels of Fig. 3. This is also clearly manifested 
in the PDFs showing decreasing significance of the distribu- 
tion of the cold phase. Due to the skewness of the density dis- 
tribution towards higher densities, also the pressure distribu- 
tion was asymmetric in the thermally stable case. In the ther- 
mally unstable cases, the pressure PDF is also clearly asym- 
metric, but in addition exhibits an extended wing towards the 
low pressures, and this wing gets more pronounced with in- 
creasing forcing amplitude, visible in the rightmost panel of 
Fig. 4. In our study, the vigor of forcing is still modest and 
the resulting Mach numbers are reasonably low in compari- 
son to, e.g. Gazol et al. (2005), who studied a similar system 
in both the sub- and supersonic regimes. In their results, the 
signatures of TI in the PDFs are much less pronounced; on the 
other hand, the extended wing seen to develop in our pressure 
PDFs, is not reported in theirs. 

As evident from Table 2, even though the cooling functions 
were scaled to produce the same integrated cooling over the 
whole temperature range for constant density, the thermal en- 
ergy content in the thermally unstable runs is roughly half of 
that of the thermally stable counterparts. Obviously, this dif- 
ference is due to the fact that the density is not constant, but 
actually very different in the thermally stable versus thermally 
unstable systems, as was discussed in length in this section, 
causing a nonlinear back-reaction to the cooling process, sig- 
nificantly enhancing its efficiency. For the thermally stable 
runs, however, the thermal energy is almost constant with any 
value of forcing investigated. For thermally unstable cases, 
the thermal energy is an increasing function of forcing. 

3.2. Non-helical forcing 

Using the same setup as for the hydrodynamic runs, we ad- 
ditionally include a very small random seed field of rms am- 
plitude Bi„it = 4x 10 _4 /iG and zero mean to the system, and 
follow the evolution of the system with the full set of MHD- 
equations with both cooling functions. In the following, we 
divide the magnetic field into its large-scale and fluctuating 
components, i.e. B tot = B + b, the condition (b)=0 hold- 
ing. In the runs TSanh-TScnh (thermally stable) and SSanh - 
SScnh (thermally unstable) the forcing function used has no 
relative helicity, i.e. (7=0. Large-scale dynamo action, i.e. 
the growth of the large-scale magnetic field B that has a non- 
zero mean, has been found in systems with non-helical forc- 
ing when large-scale shear has been included (Yousef et al. 
2008). Systems such as the one studied here, however, can 
be expected to develop only small-scale dynamo action, i.e. 
the growth of the fluctuating component b, when subjected to 
non-helical forcing. The results for the six non-helical runs 
produced are listed in Table 2. The forcing amplitude was in- 
creased as high as was permissible by the numerical scheme 
without invoking shock-capturing viscosities. No small-scale 
dynamo action was observed, i.e. the total magnetic field 
strength was exponentially decaying identically to the helical 
cases without large-scale dynamo action (see Fig. 7), in any of 
the runs, i.e. in the range of Reynolds numbers 31 — 97. This 
indicates that in the system under investigation, the critical 
Reynolds number for small-scale dynamo action is larger than 
the ones reached in our runs, and that if magnetic fields are 
generated in the helical cases investigated, they arise solely 
due to the action of a large-scale dynamo. 

Haugen et al. (2004) studied the onset of dynamo action 
in an isothermal system with non-helical forcing. As in this 
study, the vigor of forcing was increased, and it was observed 



that the critical Reynolds number for the onset of the small- 
scale dynamo was increasing with the Mach number. Shock- 
capturing viscosities were used to reach the trans- and su- 
personic regimes. For subsonic flows, the critical Reynolds 
number for small-scale dynamo action was of the order of 
35, while for supersonic flows with Mach numbers exceed- 
ing unity, the critical Reynolds number was almost doubled 
up to 70. Our weakest forcing cases with modest Mach num- 
bers have Reynolds numbers slightly below (TSanh, 31) or 
slightly exceeding the limit (SSanh, 39) found by Haugen 
et al. (2004), and in the light of their results, it is under- 
standable that no small-scale dynamo action can be seen. The 
modest and high amplitude forcings used here do not yet pro- 
duce Mach numbers even close to the supersonic limit, but 
their Reynolds numbers are clearly larger than the critical lim- 
its found by Haugen et al. (2004). Still, no small-scale dy- 
namo action is seen, indicating that the inclusion of thermo- 
dynamics with the quite complex heating and cooling proper- 
ties further alters the Mach-number dependence of the critical 
Reynolds number for small-scale dynamo action. 

As can be seen from Table 2, the non-helical runs produce 
somewhat weaker turbulence than the runs with purely heli- 
cal forcing er=l. This is also reflected by the density distribu- 
tion, so that the density contrast is reduced in the non-helical 
runs compared to the helical counterparts. This is manifested 
by smaller rms values for velocity, and thereby also smaller 
Reynolds numbers. The differences in the Reynolds numbers 
are not very large, so we believe that our conclusion about the 
absence of small-scale dynamo in the helical runs still holds. 

3.3. Thermally stable magnetohydrodynamic runs 

Next we repeat the magnetohydrodynamic runs with helical 
forcing, starting with the thermally stable cooling function, 
which is expected to produce large-scale dynamo action (e.g. 
Brandenburg 2001). The Run TSf with the forcing amplitude 
/ = 7 is sub-critical to dynamo action, i.e. the small seed 
field exponentially decays. The Run TSe with / = 10, how- 
ever, produces a very slowly growing magnetic field, with the 
time scale of growth of roughly 500Myrs. The growth of the 
field saturates only after 4 Gyrs, the saturation energy being 
less than half of the kinetic energy of turbulence. Inspecting 
the Reynolds numbers realized in Runs TSe and TSf from Ta- 
ble 2, the critical Reynolds number for the large-scale dynamo 
in this system, therefore, is roughly between 12-16. 

Runs with higher forcing (Runs TSa, TSd, TSb and TSe), 
with a forcing ranging between 20-70, all produce dynamo 
action that leads to the generation of magnetic fields, the en- 
ergy of which is in rough equipartition with the kinetic en- 
ergy of turbulence for intermediate forcings (TSa and TSd), 
but slightly lower than the equipartition value for the high- 
est forcings investigated (TSb and TSe). The growth rate 
and saturation strength of the field are increasing as func- 
tion of increasing forcing, and thereby also as functions of the 
Reynolds number Rm, so that for the highest forcing investi- 
gated, the time scale of growth is roughly 27 Myrs, building 
up the magnetic field in a few hundreds of Myrs. A simi- 
lar, but isothermal, system has been investigated earlier by 
Brandenburg (2001), who found otherwise similar behavior, 
but the energy of the magnetic fields generated in the sys- 
tem exceeded the kinetic energy of turbulence. The build-up 
of the large-scale field from the equipartition to the super- 
equipartition values was observed to occur on the diffusive 
timescale, and therefore be dependent on Rm. Such behavior 
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Model f h ^imz gig; ^i~ y ^(^1 vm)Z 

TSe 10 1 - O06 O05 007 O05 - 

TSa 20 1 - - 0.49 - 0.43 

TSd 35 1 0.70 .... 0.79 

TSb 50 1 0.98 .... 1.01 

TSc 70 1 - L18 - L23 - - 

SSf 15 1 - 0.09 - 0.09 

SSa 20 1 - - 0.21 - 0.23 

SSb 50 1 0.58 .... 0.60 

SSc 70 1 - 0.80 - 0.81 

TABLE 3 

The rms-stregth of the k = 1-mode, measured in ^G, generated in the runs. 
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FIG. 6. — The calculated flux diagnostic from Eq.(12) for the thermally 
stable run TSa (dashed-dotted), TSd (dashed) and TSc (solid line). 

is not seen in the system of this study. 

Next, we study the spatial distribution of the magnetic 
fields, and convince ourselves that the dynamo seen in the 
helical forcing case is indeed of the large-scale type, being 
capable of generating magnetic fields on spatial scales larger 
than the forcing scale of the turbulence, kf/k\ = 3. In the 
two leftmost panels of Fig. 5, we plot two components of the 
magnetic field for Run TSa after the growth of the magnetic 
field has saturated. For the B y component, a clear sinusoidal 
wave of the Fourier mode k = 1 over the x-coordinate is seen, 
while B z shows the same mode with a cosinus dependence of 
the spatial coordinate. The patterns seen are very similar to 
the study of Brandenburg (2001), and the field is clearly co- 
herent over larger scales as the forcing scale. Such a mean 
field vanishes when averaging over the whole volume, but the 

rms-strength of the large-scale component \J {By) yz can be 

recovered if averaging over y- and z-directions, denoted with 
() , is used; we use averaging of this type in Table 3 to esti- 
mate the generated large-scale field strengths. 

The scales of the large-scale fields are the same for all the 
runs independent of forcing, but the direction of the gener- 
ated mean field varies almost from one run to another. As can 
been seen from Table 3, the weakest forcing case, Run TSf, 
the Reynolds number being very close to the marginal one 
and the time scale of growth very long, exhibits mean fields 
in every direction. Run TSa, the magnetic field shown in 
Fig. 5, develops strong mean components of B y and B z in 
the x-direction. The intermediate forcing cases, Runs TSd 



and TSb, on the other hand, show mean components of B x 
and B z over the y-direction. Finally, the strongest forcing 
case, Run TSc, develops the strongest mean fields in B x and 
By over the z-direction. Also in the isothermal study of Bran- 
denburg (2001), any direction of the mean field could be pre- 
ferred, only depending on the fine details of the initial random 
distribution. As is evident when comparing the rms values of 
the total magnetic field, presented in Table 2, to the strength 
of the mean component in Table 3, almost all magnetic energy 
is contained in the k = 1 mode. As the forcing is increased, 
small contribution of the small-scale field is generated, but in 
all cases investigated this contribution is negligible, verifying 
the assumption that no small-scale dynamo is present in the 
models. 

When calculating the flux diagnostic, defined by Eq. (12), 
for a magnetic field with a mean component at the wavenum- 
ber k = 1 over a certain direction, having different signs 
and roughly equal magnitudes in the different halves of the 
box, the total magnetic flux would evidently be close to zero. 
Therefore, we need to use absolute values of the field strength, 
\B\, in Eq. 12. Furthermore, as the distribution of the mean 
fields can be quite complex, we calculate an average value (if?) 
of <f? xy , <f? yz , and f IZ . In Fig. 6 we show plots of the mag- 
netic flux diagnostic for the thermally stable runs developing 
dynamo action. For the weakest forcing case, Run TSa, the 
magnetic and density fields show signs of anti-correlation for 
low densities, n < lcm -3 , and no clear correlation for den- 
sities above that limit. The anti-correlation at lower densities 
is observed to decrease with increasing forcing and dynamo 
efficiency, but only for the highest forcing and most efficient 
dynamo, Run TSc, are signs of positive correlation at higher 
densities, 10-35 cm -3 , established. This clearly shows that 
the compression due to the forcing alone is very inefficient 
in producing any significant positive correlation between the 
density and magnetic field. 

3.4. Thermally unstable magnetohydrodynamic runs 

Finally, we produce magnetohydrodynamic runs with heli- 
cal forcing using the thermally unstable cooling function, the 
results being listed in the lower part of Tables 2 and 3 and the 
time evolution of the magnetic energy for some runs shown in 
Fig. 7. As is evident from this table, when comparing these 
runs with their thermally stable counterparts with the same 
magnitude of forcing, it can be noted that the turbulence is 
more vigorous with the thermally unstable cooling function 
in the sense that rms values of velocity and Reynolds num- 
bers are larger for any fixed forcing amplitude. The kinetic 
energy, however, is smaller for the thermally unstable runs 
than for the thermally stable ones, explained by the density 



9 



1 o 




10"' 



0.0 0.2 0.4 0.6 0.8 1.0 

t [Gyr] 

FIG. 7. — Time evolution of the magnetic fields in some of the thermally 
unstable runs. 
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FIG. 8. — Left panel: Thermal energies and the ratio of magnetic to kinetic 
energy as function of magnetic Reynolds number for runs with different cool- 
ing functions. Stars: thermal energy in the thermally stable MHD Runs TSf 
- TSc. Crosses: thermal energy in the thermally unstable MHD Runs SSf - 
SSc. Diamonds: magnetic to kinetic energy ratio in Runs TSf— TSc. Trian- 
gles: magnetic to kinetic energy ratio in Runs SSf— SSc. Right panel: Den- 
sity contrast Ap from MHD runs (stable with stars, unstable with diamonds) 
compared with their hydrodynamic counterparts (stable with crosses, unsta- 
ble with traingles) as function of the forcing amplitude, /. 



distribution being dominated by the warm diffuse gas with 
low density for the thermally unstable setup. 

For the Run SSf with the forcing amplitude /=15 producing 
turbulence with Rm « 45, a very slowly growing dynamo is 
found, as seen from Fig. 7, whereas the magnetic field of the 
Run SSe with Rm « 28 is seen to decay exponentially. The 
critical Reynolds number, therefore, lies somewhere in be- 
tween 28^-5, roughly twice the value found for the thermally 
stable runs. In other words, the large-scale dynamo is harder 
to excite in the thermally unstable system. As is evident from 
the left panel of Fig. 8, where we plot the dependence of ther- 
mal energy and the ratio of magnetic to kinetic energy ratio on 
the magnetic Reynolds number Rm, the magnetic energy in 
the saturated state remains below the kinetic energy of turbu- 
lence up to the forcing amplitude /=35 in Run SSd, for which 
the Reynolds number is 58. For the thermally stable runs the 
dynamo efficiency of this level was reached with /=20 with 
Rm=27, the difference being again by a factor of two. For the 



strongest forcings investigated (Runs SSb and SSc), the mag- 
netic energy exceeds the kinetic energy of turbulence. Such 
behavior was not seen in the thermally stable runs. 

The density contrast can be seen to become affected by the 
presence of the magnetic field; this effect is visible in Fig. 8 
right panel, where we plot Ap both for the MHD runs and 
their hydrodynamic counterparts. For both the stable and un- 
stable runs the trend appears very similar: as long as the mag- 
netic field energy remains clearly below the kinetic energy of 
turbulence, i.e. the dynamo efficiency is still low, the density 
contrast becomes enhanced. For forcings capable of generat- 
ing a magnetic field close to or exceeding the kinetic energy of 
turbulence (Runs TSc, SSb, and SSc) the density contrast be- 
comes reduced. This can be interpreted as the lowered capa- 
bility of either the turbulence or the condensation mode of the 
thermal instability to create density structures in the presence 
of a strong dynamo-generated magnetic field, in agreement 
with the results from earlier magnetohydrodynamic studies 
without self-consistent dynamo action, e.g. by Ostriker et al. 
(1999) and Padoan & Nordlund (1999). 

In Fig. 9 we show PDFs of the magnetohydrodynamic runs 
with the strongest forcings (Runs SSb and SSc) compared to 
the hydrodynamic counterparts (Runs HDSSb and HDSSC). 
The main effect due to the presence of the magnetic field seen 
in them is the decreased amount of dense gas at the very high 
end of the density distribution, and a slight increase of very 
low density gas for Run SSc, while hardly any difference can 
be seen in the temperature distributions. The pressure PDF, 
therefore, reflects the decreased densities, the higher pressure 
wing being reduced for both Runs SSb and SSc. This implies 
that the stronger the magnetic field grows, the more it resists 
the cloud compression due to thermal instability. In a similar 
study, Gazol et al. (2009) reported that an initially uniform 
magnetic field of varying strengths caused an extended low- 
pressure wing to their pressure distribution; in our case, the 
wing was seen already in the hydrodynamic regime, and the 
effect of the magnetic field on it is barely visible. In a re- 
lated, but not completely comparable, study by de Avillez & 
Breitschwerdt (2005), the effect of the presence of a magnetic 
field on a supernova-regulated flow was studied. In contrast 
to our results, the amount of cold and dense gas was observed 
to increase in the magnetohydrodynamic versus the hydrody- 
namic setup. The differences seen between our model to the 
other two might be connected to the fact that in our case self- 
consistently dynamo-active flows were studied in contrast to 
imposed ones passively advected by the flow. 

Finally, we plot the flux diagnostic (ty) for the dynamo- 
active thermally unstable runs in Fig. 10. If compared to the 
corresponding plot for thermally stable runs, Fig. 6 showing 
uncorrected density and magnetic fields except for the high- 
density end of the run with the highest forcing, the action 
of thermal instability on the system produces clear positive 
correlation between the two quantities. The correlation for 
all the cases investigated is considerably weaker, of the form 
B oc p° 2 (thick line plotted in the figure), compared to the 
one observed in the interstellar matter. The strongest correla- 
tion is seen for the Run SSd, the trend declining somewhat for 
higher forcings. For all the runs, there seems to be a tendency 
for the correlation to break down at the highest densities. 

4. CONCLUSIONS 

In this paper we have investigated the role of thermal in- 
stability for large-scale dynamo action in galaxies and in pro- 
ducing the density-magnetic field correlation in the interstel- 
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FIG. 9. — Probability density functions for the thermally unstable hydrodynamic runs HDSSc (solid line) compared with the magnetohydrodynamic counterpart 
SSc (dashed line), and HDSSb (dotted line) with SSb (dashed-dotted line). On the left we show the PDF of logarithmic density, in the middle the PDF of 
logarithmic temperature, and on the right the PDF of logarithmic pressure. 

turbulent velocity field. In terms of the rms-magnetic field 
strength, the thermally stable runs produced roughly 1.5 times 
larger magnetic fields than the thermally unstable counter- 
parts per fixed forcing, the maximal values being of the order 
of 1.5/iG for the strongest forcing. The dynamos, no matter 
which cooling function was used, always produced magnetic 
fields with a mean component at the wavenumber k=l. The 
large-scale field could be generated in any direction, and show 
systematic variation over any spatial coordinate, depending on 
the strength of the forcing. 

In the thermally stable runs the denser structures were 
formed due to the action of the turbulent velocity field itself. 
The enhanced density contrast produced by increasing forc- 
ing, however, was found to be very inefficient in producing 
a density-magnetic field correlation in the system. In fact, 
a positive correlation was found only for the highest forcing 
case investigated, and even then only at the limit of the high- 
est densities. For all the other cases, our flux diagnostic in- 
dicated a tendency for uncorrelation between the quantities. 
For the thermally unstable dynamo-active cases investigated, 
a positive correlation of the form B oc p° 2 was observed to 
be established, but a tendency of this rather weak correlation 
to break down at the high density limit was seen. It can be 
concluded that an instability of any 'condensation' type, i.e. 
the gravitational or condensation mode of thermal instability, 
is instrumental in establishing density-magnetic field correla- 
tion in the ISM, but based on this study, the correlation pre- 
dicted due to thermal instability is weak, and especially so for 
high densities. 
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FIG. 10. — The calculated flux diagnostic from Eq.(12) for the ther- 
mally unstable runs SSd (dashed-dotted), SSb (dashed) and SSc (solid line). 
The thickest line shows a fit of the type B oc p 2 . 

lar matter. A periodic cubic magnetohydrodynamic setup with 
external helical forcing was used to establish large-scale dy- 
namo action in the system, which was subject either to ther- 
mally stable or unstable cooling functions having the same 
integrated net cooling for constant density. Reference runs 
without magnetic fields and with non-helical forcing were 
made, the former to isolate the effect of the magnetic fields 
on the system, and the latter for detecting if small-scale dy- 
namo action was produced in the Reynolds number regime 
investigated. The latter set of runs showed that no small-scale 
dynamo action was to be expected in the helical runs, the crit- 
ical Reynolds number for small-scale dynamo action in this 
particular system being larger than 97. This number is larger 
than obtained in similar, but isothermal, setups investigated. 

In the thermally stable runs, the critical Reynolds num- 
ber of large-scale dynamo action was calculated to be in be- 
tween 12-16, whereas for the thermally unstable setup the 
corresponding range was found to be roughly twice that, 28- 
45; the thermal instability, therefore, makes the onset of the 
large-scale dynamo harder. In contrast to the thermally stable 
cases investigated, that were observed to produce magnetic 
fields in equipartition with the kinetic energy of turbulence, 
the thermally unstable ones were capable of producing mag- 
netic fields with somewhat larger energy than contained in the 
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